Persistence exponents for fluctuating interfaces 
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Numerical and analytic results for the exponent 6 describing the decay of the first return proba- 
bility of an interface to its initial height are obtained for a large class of linear Langevin equations. 
The models are parametrized by the dynamic roughness exponent 0, with < < 1; for = 1/2 
the time evolution is Markovian. Using simulations of solid-on-solid models, of the discretized 
continuum equations as well as of the associated zero-dimensional stationary Gaussian process, we 
address two problems: The return of an initially flat interface, and the return to an initial state with 
fully developed steady state roughness. The two problems are shown to be governed by different 
exponents. For the steady state case we point out the equivalence to fractional Brownian motion, 
which has a return exponent 9s = 1 — 0- The exponent 6q for the flat initial condition appears to 
be nontrivial. We prove that 6 -> oo for -» 0, 6 > Os for < 1/2 and 6 < 6 S for > 1/2, 
and calculate #o,s perturbatively to first order in an expansion around the Markovian case = 1/2. 
Using the exact result 8s ~ 1 — 0, accurate upper and lower bounds on 8o can be derived which 
show, in particular, that do > (1 — 0) 2 /0 for small 0. 

PACS numbers: 02.50.-r, 05.40.+j, 81.10.Aj 



I. INTRODUCTION 

The statistics of first passage events for non-Markovian 
stochastic processes has attracted considerable recent in- 
terest in the physical literature. Such problems appear 
naturally in spatially extended nonequilibrium systems, 
where the dynamics at a given point in space becomes 
non-Markovian due to the coupling to the neighbours. 
The asymptotic decay of first passage probabilities turns 
out to be hard to compute even for very simple systems 
such as the one-dimensional Glauber model jjj or the 
linear diffusion equation with random initial conditions 
Indeed, determining the first passage probability of 
a general Gaussian process with known autocorrelation 
function is a classic unsolved problem in probability the- 
ory §§. 

In this paper we address the first passage statistics of 
fluctuating interfaces. The large scale behaviour of the 
models of interest is described by the linear Langevin 
equation 



A(t, t') ee (h(x, t)h(x, t')} = K[(t' + t) 2f3 - \t' - t\ 



|2/3l 



(2) 



dh 



-(-V 2 ) z/2 h + v 



(1) 



for the height field h(x,i). Here the dynamic exponent 
z (usually z = 2 or 4) characterizes the relaxation mech- 
anism, while rj{x,t) is a Gaussian noise term, possibly 
with spatial correlations. We will generally assume a flat 
initial interface, h(x, 0) = 0. Since (0) is linear, h(x, t) is 
Gaussian and its temporal statistics at an arbitrary fixed 
point in space is fully specified by the autocorrelation 
function computed from (Q), 



where K is some positive constant, and f3 denotes the 
dynamic roughness exponent, which depends on z and 
on the type of noise considered. For example, for un- 
correlated white noise f3 = (1/2) [1 — d/z] for a d- 
dimensional interface, while for volume conserving noise 
(3 = (1/2) [1 - (d + 2)/z] [§. An interface is rough if 
[3 > 0. In the present work we regard j3 as a continuous 
parameter in the interval ]0, 1[. Note that for j3 = 1/2 
eq.(|^) reduces to the autocorrelation function of a ran- 
dom walk, corresponding to the limit z — > oo (no relax- 
ation) of eq.(|l|) with uncorrelated white noise. 

To define the first passage problems of interest, con- 
sider the quantity 

P(t ,t) = Prob[h(x, s) ^ h(x,t ) Vs : t < s < t + t]. 

(3) 

We focus on two limiting cases. For to = 0, P(to,t) re- 
duces to the probability po(t) that the interface has not 
returned to its initial height h = at time t. This will be 
referred to as the transient persistence probability, char- 
acterized by the exponent (9q, 



(4) 



On the other hand, for to — > oo the interface develops 
roughness on all scales and the memory of the flat initial 
condition is lost. In this limit P(tQ, t) describes the return 
to a rough initial configuration drawn from the steady 
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state distribution of the process, and the corresponding 
steady state persistence probability ps(t) decays with a 
distinct exponent 9s, 



p s (t) = lim P(t ,t) 

t()—>OC 



(5) 



In general, one expects that P(to,t) ~ t~ 6s for t <C to 
and P(t ,t) ~ t~ e ° for t ^> t , with a crossover function 
connecting the two regimes. 

A particular case of the steady state persistence prob- 
lem was studied previously in the context of tracer diffu- 
sion on surfaces [ ?J . In this work it was observed that the 
distribution of first return times has a natural interpre- 
tation as a distribution of trapping times during which 
a diffusing particle is buried and cannot move; thus the 
first passage exponent may translate into an anomalous 
diffusion law. A simple scaling argument (to be recalled 
below in Section V) was used to derive the relation 



1-0, 



(6) 



which was well supported by numerical simulations for 
= 1/8. A primary motivation of the present work 
is therefore to investigate the validity of this relation 
through simulations for other values of and refined an- 
alytic considerations, as well as to understand why it fails 
for the transient persistence exponent 9o- 

The paper is organized as follows. In the next section 
we convert the non-stationary stochastic process h(x,t) 
into a stationary Gaussian process in logarithmic time 
This representation will provide us with a number 
of bounds and scaling relations, and will be used in the 
simulations of Section IV. C. A perturbative calculation 
of the persistence exponents in the vicinity of = 1/2 is 
presented in Section III. The simulation results are sum- 
marized in Section IV. Section V reviews the analytic ba- 
sis of the relation (||) and makes contact to earlier work 
on the return statistics of fractional Brownian motion, 
while Section VI employs the expression (||) for #5 to nu- 
merically generate exact upper and lower bounds on 9q. 
Finally, some conclusions are offered in Section VII. 



II. MAPPING TO A STATIONARY PROCESS 



Following Refs. 



we introduce the normalized ran- 



dom variable X = hj \J (h 2 ) which is considered a func- 
tion of the logarithmic time T = Int. The Gaus- 
sian process X(T) is then stationary by construction, 
(X(T)X(T')) = /o(T-T'), and the autocorrelation func- 
tion fo obtained from (||) is 



fo(T) = cosh(T/2) 2 ' 3 - I sinh(T/2)| 2 ' 3 . 



(7) 



In logarithmic time the power law decay (^) of the 
persistence probability becomes exponential, Po(T) ~ 
exp(— 9qT), and the task is to determine the decay rate 
60 as a functional of the correlator f . 



Similarly a normalized stationary process can be as- 
sociated with the steady state problem. First define the 
height difference variable 



H(x, t; t ) = h(x, t + to) — h(x, t ) 



(8) 



and compute its autocorrelation function in the limit 

t Q — * 00, 

A s {t,t')= lim (H(x,t;t )H(x,t';t )) 

= lim [A(t +t,t + t') - A{t + t, t ) 

to — >oo 

-A(t ,t + t')+A(t ,to)] 

= K[t 20 + t' 2 ? -\t' -t\ 2 % (9) 

which is precisely the correlator of fractional Brownian 
motion with Hurst exponent ^ (see Section V). Next 
As(t, t') is normalized by As{t, t)As(t', t') and rewrit- 
ten in terms of T — Int. This yields the autocorrelation 
function 

f s (T) = cosh(/3T) - i|2sinh(T/2)| 2 ' 3 . (10) 

Comparison of eqs.(^) and (|l^) makes it plausible that 
the two processes have different decay rates of their per- 
sistence probabilities. Both functions have the same type 
of short time singularity 



fo,s(T) = 1 ~ 0(\T 



2/3 \ 



T ^ 0, 



(11) 



which places them in the class a = 20 in the sense of 
Slepian ||. However, for large T they decay with differ- 
ent rates, /o,s(T) ~ exp(— \o,sT) for T — > 00, where 



X o = l-0, X s = mm[0,l-0] 



(12) 



can be interpreted, in analogy with phase ordering ki- 
netics, as the autocorrelation exponents || of the two 
processes. 

For a stationary Gaussian process with a general au- 
tocorrelator f(T), the calculation of the decay exponent 
of the persistence probability is very hard. Only in a 
very few cases exact results are known ||. Approximate 
results can be derived for certain classes of autocorre- 
lators f(T). For example, when f(T) = 1 - 0(T 2 ) for 
small T (an example being the linear diffusion equation 
0] ) , the density of zero crossings is finite and an indepen- 
dent interval approximation (IIA) Q gives a very good 
estimate of 9. However, for any other process for which 
f(T) = 1 - 0(\T\ a ) for small T with a < 2, the density 
of zeros is infinite and the IIA breaks down. For gen- 
eral proceses with a = 1, a perturbative method (when 
the process is not far from Markovian) and an approxi- 
mate variational method was developed recently |5|] . This 
method will be applied to the present problem in Section 
III. In the remainder of this section we collect some exact 
bounds on 8; further bounds will be derived in Section 
VI. 
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Slepian || has proved the following useful theorem for 
stationary Gaussian processes with unit variance: For 
two processes with correlators fi(T) and f2(T) such that 
fi(T) > h{T) > for all T, the corresponding persis- 
tence probabilities satisfy Pi(T) > p2(T); in particular, 
the inequality 9\ < 02 holds for the asymptotic decay 
rates. By applying this result to the correlators (0) and 
( |l0| ) we can generate a number of relations involving the 
return exponents O and 9s- For example, taking the 
derivative of (Q) with respect to (3 one discovers that 
/ (T) increases monotonically with decreasing (3 for all 
T, and consequently 

0o((3) > ((3') if P<&. (13) 
For (3 < (3 1 < (21n2) _1 the stronger inequality 

(l-(3')0o(f3)>(l-P)6o(P') (14) 

is proved in the Appendix. 

Moreover, rewriting (fiXI) in the form 



f s (T) = e-^ + ie^ r l[l - e- 2 ^ T l - (1 - e -l T l) 2 ^] 



(15) 

it is evident that f s (T) < exp(-/3|T|) for (3 < 1/2 and 
f s (T) > exp(-/3\T\) for (3 > 1/2. A process charac- 
terized by a purely exponential autocorrelation function 
exp(— A|T|) is Markovian, and its persistence probability 
can be computed explicitly ; the asymptotic decay rate 
is equal to the decay rate A of the correlation function. 
Thus the fact that fs(T) can be bounded by Markovian 
(exponential) correlation functions supplies us with the 
inequalities 



0s > A 13 < 1/2 
0s<f3, 13 > 1/2. 
The last inequality can be sharpened to 



<l/2, /3 > 1/2. 



(16) 



(17) 



This will be demonstrated in the Appendix, where we 
also prove that 



and 



#0 > 1-/3 for (3 < 1/2 
%< 1-/9 for (3> 1/2 

6> > S for (3 < 1/2 
O < S for /3 > 1/2. 



(18) 



(19) 



Next we record some relations for special values of (3. 
We noted already that for (3 — 1/2 the interface fluc- 
tuations reduce to a random walk, corresponding to the 



Markovian correlator fo(T) — fs(T) = 
which 6 = A = 1/2 g. Hence 

o (l/2) = 05(1/2) = 1/2. 



exp(-|T|/2), for 



(20) 



For (3 = 1 both (g) and (|l0|) become constants, f Q (T) = 
fs(T) = 1. This implies that the corresponding Gaussian 
process is time-independent, and consequently 



lim O 



lim 9 s (f3) = 0. 



(21) 



For (3 — > the transient correlator (|7|) degenerates to 
the discontinuous function /o(0) = 1, ^(T 1 > 0) = 0. 
Since this is bounded from above by the Markovian cor- 
relator f(T) = e _A ' T l for any A, we conclude that 



lim O 09) 



00. 



(22) 



In contrast, the steady state correlator tends to a nonzero 
constant, lim^o fs(T) = 1/2 for T > 0, with a discon- 
tinuity at T = 0, and therefore 9s is expected to remain 
finite for (3 — ► 0. Note that all the relations derived for (9$ 
- equations ( ^6| , ^7 21) - are consistent with 63 = 1 — 0. 



III. PERTURBATION THEORY NEAR = 1/2 

We have already remarked that both the steady state 
and the transient processes reduce to a Markov process 
when (3 — 1/2. Two of us have developed a perturba- 
tion theory for the persistence exponent of a stationary 
Gaussian process whose correlation function is close to 
a Markov process ||. When the persistence probability 
for a Markov process is written in the form of a path 
integral, it is found to be related to the partition func- 
tion of a quantum harmonic oscillator with a hard wall 
at the origin. The persistence probability for a process 
whose correlation function differs perturbatively from the 
Markov process, i.e. whose autocorrelation function is 



/(T)=exp {-\\T\) + ecj>{T), 



(23) 



may then be calculated from a knowledge of the eigen- 
states of the quantum harmonic oscillator. In ( |23| ) we 
have used the same normalization, /(0) = 1, as elsewhere 
in this paper. With this normalization, (f>(0) = 0. (Note 
that a different normalization was employed in ref. ||.) 

If [3 = 1/2 + e, equations (|7|) and ( |l0| ) may be written 
in the form 



1 \T\ 
fo,s = cxp ( — 



■e<Ks(|T|) + 0(e 2 ), 



(24) 



where 



T ( T 
2 cosh — In cosh — 
2 V 2 



T ( 

2 sinh - In f sinh - ) , (25) 



4>s 



T ( ( T 

sinh — T - 2 In 2 sinh - 
2 V V 2 



(26) 
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The result for the persistence exponent (equivalent to 
equation (7) of reference ||) may most conveniently be 
written in the form 

9 = \{l-e^f (f)(T) [1 - exp(-2AT)r 3/2 dT J . 

(27) 

Substituting ( p5| ) and (|2^) into (^), one hnds (after some 
algebra) 

9 = ^~e(2V2-l) + O(e 2 ), (28) 
e s = 1 --e + 0{e% (29) 

Eqn. ( p9| ) agrees with the relation #s = 1 — /3, while 
( p8| ) compares favourably with the stationary Gaussian 
process simulations for 8 = 0.45 and j3 = 0.55 to be 
presented in Section IV. C. 



IV. SIMULATION RESULTS 

A. Solid-on-solid models 

Simulations of one-dimensional, discrete solid-on-solid 
models were carried out for 8 = 1/8, /3 = 1/4 and 
8 = 3/8. The case (3 = 1/8 describes an equilibrium 
surface which relaxes through surface diffusion, corre- 
sponding to z — 4 in (Q) and volume conserving noise 
with correlator 

(r)(x, t)rj(x', t')) = -V 2 5(x - x')5{t - t'). (30) 

The cases (3 — 1/4 and 8 = 3/8 are realized for noncon- 
served white noise in ( |lh and dynamic exponents z = 2 
and z = 4, respectively || . 

In all models the interface configuration is described 
by a set of integer height variables hi defined on a one- 
dimensional lattice i = 1,...,L with periodic boundary 
conditions. For simulations of the transient return prob- 
lem, large lattices [L = 2 x 10 5 - 2 x 10 6 ) were used, 
while for the steady state problem we chose small sizes, 
L = 100 — 200 for dynamic exponent z — 4 and L = 1000 
for z — 2, in order to be able to reach the steady state 
within the simulation time. 

The precise simulation procedure is somewhat depen- 
dent on whether the volume enclosed by the interface is 
conserved (as for (3 = 1/8) or not. To simulate the tran- 
sient return problem with conserved dynamics, the inter- 
face was prepared in the initial state hi = 0, and each site 
was equipped with a counter that recorded whether the 
height hi had returned to hi = 0. The fraction of coun- 
ters still in their initial state then gives the persistence 
probability po(t)- For the steady state problem the inter- 
face was first equilibrated for a time t cq large compared 



to the relaxation time ~ L z HJ!^]. Then the configu- 
ration hi(t C q) was saved, and the fraction of sites which 
had not yet returned to hi(t eq ) was recorded over a pre- 
scribed time interval t cq < t < £ e q +£i - At the end of that 
interval the current configuration hi(t cq + ti) was chosen 
as the new initial condition, and the procedure was re- 
peated. After a suitable number of repetitions (typically 
10 4 for z = 4, and 2000 for z = 2), the surviving fraction 
gives an estimate of ps- 

The models used in the cases 8=1/4 and 8 = 3/8 are 
growth models, in which an elementary step consists in 
chosing at random a site i and then placing a new par- 
ticle, hj — > hj + 1, either at j = i or at one of the two 
nearest neighbour sites j = i ± 1, depending on the local 
environment. For these nonconserved models the pro- 
cedures described above have to be modified such that 
the calculation of the surviving fractions po and p$ is 
performed only when a whole monolayer - that is, one 
particle per site - has been deposited. At these instances 
the average height is an integer which can be subtracted 
from the whole configuration in order to decide whether 
a given height variable has returned to its initial state 
when viewed in a frame moving with the average growth 
rate. 

We now briefly describe the results obtained in the con- 
served case. In JjJ the steady state return problem for 
(3 =1/8 was investigated in the framework of the stan- 
dard one-dimensional solid-on-solid model with Hamilto- 
nian 

H = jJ2\hi+i-hi\ (31) 



and Arrhenius-type surface diffusion dynamics We 
have extended these simulations to longer times and to 
different values of the coupling constant J. Figure 1 
shows that the exponent 9$ is independent of J, and that 
its value is numerically indistinguishable from 9g = 7/8 
predicted by 

Since the transient persistence probability po decays 
very rapidly for 8 — 1/8, a more efficient model was 
needed in order to obtain reasonable statistics. We there- 
fore used a restricted solid-on-solid (RSOS) model in- 
troduced by Racz et al. E^). In this model the nearest 
neighbour height differences are restricted to 

\h i+1 -h i \<2. (32) 

In one simulation step a site i is chosen at random, and 
a diffusion move to a randomly chosen neighbour is at- 
tempted. If the attempt fails due to the condition (|32"|), a 
new random site is picked. Figure 2 shows the transient 
persistence probability po(t) obtained from a large scale 
simulation of this model. The curve still shows consid- 
erable curvature, and we are only able to conclude that 
probably 9 > 3.3 for this process. 
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FIG. 1. Steady state persistence probability ps(t) for 
the Arrhenius surface diffusion model with coupling constant 
J = 0.25, 0.5 and 1. Systems of size L = 100 were equili- 
brated for t cq = 5 x 10 7 attempted moves per site. Then data 
were collected over 10 4 time intervals of length t\ = 10 4 . The 
dashed line has slope 6s = 1- /3 = 7/8. 



Figure 2 also shows transient results for j3 = 1/4 and 
(3 = 3/8. In the former case we used a growth model in- 
troduced by Family JL5) , in which the deposited particle 
is always placed at the lowest among the chosen site i 
and its neighbours, whereas for (3 — 3/8 we used the cur- 
vature model introduced in Jl6| . Our best estimates of Oq 
for these models are collected in Table I, along with the 
values for 9 s which agree, within numerical uncertainties, 
with the relation m) in all cases. 



B. Discretized Langevin equations 



We solved equation (|lj) in discretized time and space 
for the real valued function h(xi,t n ), where t n = nAt 
and Xi = iAx with n = 0, 1, 2, . . . and i = 0, . . . , L — 1 in 
a system with periodic boundary conditions. 

For the time discretization we used a simple forward 
Euler differencing scheme |l7j : 

dh(x u t n ) _ h(xj,t n+ i) - h{xj,t n ) 

dt ~ At ' [ ' 
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The spatial derivatives for the cases z — 2 and z = 4 
considered in the simulations were discretized as 



V 2 h{xi) = h^i) ~ 2h( Xi ) + h(x l+1 ), (34) 
(V 2 ) 2 ^) ee h( Xl - 2 ) ~ kh{xi-i) + 6h(xi) 

-4h(x i+1 ) + h(x i+2 ) (35) 



for the function h(xi) = h(xi,t n ) at any given time t n . 
Here and in the simulations, the spatial lattice constant 
Air was set to unity. 

With these definitions, we iterated the equation 



h(xi,t n +i) = h{x h t n ) - At(~S7 2 y/ 2 h{ Xl , t n ) 

+ aVAtr){xi,t n ) (36) 



FIG. 2. Transient persistence probability po(t) for three 
different solid-on-solid models described in the text. The sys- 
tem sizes used were L — 2 x 10 5 for j3 — 1/8, and L = 2 x 10 6 
for (3 = 1/4 and (3 = 3/8. For /3 = 1/8 (/3 = 1/4) an average 
over 1000 (10) runs was taken, while the data for /3 = 3/8 
constitute a single run. The slopes of the dashed lines cor- 
respond to the exponent estimates in the fourth column of 
Table I. 



where rj(xi,t n ) is a Gaussian distributed random number 
with zero mean and unit variance whose correlations will 
be specified below. 

Von Neumann stability analysis (lj] shows that values 
At < 1/2 and 1/8 for z — 2 and 4, respectively, have 
to be used to keep the noise-free iteration stable. The 
simulations showed that the scheme remained stable with 
At = 0.4 and 0.1 even in the noisy case. 
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FIG. 3. (a) Steady state persistence probability ps(t) from 
the numerical solution of the discretized Langevin equation 
with 2 = 2 and correlated noise (p = —0.1,0.1,0.3). The 
thick lines represent fits to the last decade of ps{t). The 
slopes are given in the second column of Table I. (b) Same as 
(a) for the transient persistence probability po(t). 

We used white noise rj w (xi, t n ) with a correlator 

{Tj w (Xi,t n )Tj w ( y Xj,t ra )) — 3i.j3n,m (^'O 



r] p (k,t) = J\S p (k)\r k exp(2m<pk) 



(40) 




with r/;. being a Gaussian distributed amplitude with zero 
mean and tpk G [0, 1[ a uniformly distributed random 
phase |l8). Due to the regularization g(0) = 1, which 
fixes the average value of S p , one has to use the modulus 
To| ) as S p can be negative for some k. Iterating eqn. 
3|) with z = 2, the correlated noise ( fh5| ) with p < 1/2 
leads to surface roughness with a measured roughness ex- 
ponent (3 m that agrees with the prediction (3 — (l + 2p)/4 
of the continuum equation (jl]) within error bars. The case 
[3 > 1/2, i.e. p > 1/2 is not accessible by this method. 

The simulated systems had a size of L = 4096, noise 
strength a — 1 and averages were typically taken over 
3000 independent runs. 

In all cases, the simulation was started with a flat inital 
condition h(x, 0) = 0. To measure the persistence prob- 
abilities, the configuration h(x,to) and the consecutive 
one h(x, to + At) were kept in memory during the simu- 
lation. In each following iteration t n = (to + At) + nAt, 
n = 1,2,3,..., an initially zeroed counter at each site 
x was increased as long as sgn(h(x,t n ) — h(x,to)) — 
sgn(h(x,to + At) — h(x,to)). The fraction of counters 
with a value larger than t then gave the persistence prob- 
ability p(t). For measurement of 9q, to was chosen to be 
zero, for 9s, At <C to <C L z , and the power law behaviour 
in the regime At <c t <c t was used. 

For comparison, 9s was also measured in small sys- 
tems L — 128 in the steady state t ^S> L z for z = 2 with 
uncorrelated noise. The results agreed with the measure- 
ments for At <§; to <C L z . However, in the steady state, 
one has to take care of measuring the power law decay 
of ps(t) only up to the correlation time ~ L z , so that 
we preferred measurement in the regime At <C <o "C L z 
here. 

Figure 3 shows typical curves ps(t) and po(t) obtained 
from the numerical solution of the discretized Langevin 
equation ( p6| ) with correlated noise for the values p = 
—0.1, 0.1 and 0.3. A summary of all measured persis- 
tence exponents 9o an d 9s as a function of the roughness 
exponent (3 can be found in Table I. 



C. Simulation of the stationary Gaussian process 



in the simulations, as well as spatially correlated noise 
r} p {xi,t n ) with 



(rip(xi,t n )rip(xj,t m )) = <J 2 g p (xi - x 3 )5 r , 



where 



Qp\X% Xj) 



1 2/>-l 



1 ,i = j 



(38) 



(39) 



and p < 1/2 a real number. A different choice of reg- 
ularizing g(0) jl8| did not change the results. If S p (k) 
denotes the discrete Fourier transform of g p , we defined 



Since a Gaussian process is completely specified by its 
correlation function, it is possible to simulate it by con- 
structing a time series that possesses the same correlation 
function. This is most easily performed in the frequency 
domain. 

Suppose fftuj) is (the Fourier transform of) a Gaussian 
white noise, with (fj(ui)fj(u>')) — 2ttS(uj + lu'). Let 



X(u)=rj(u)<Jf(u), 



(41) 



where f(uj) is the Fourier transform of the desired corre- 
lation function (notice that f(u>) is the power spectrum 



G 



of the process, so it must be positive for all u>) 
correlation function of X is 



{X(cj)X(w')) = 2nf(Lo)8{uj + u)'). 



Then the 



(42) 



The inverse Fourier transform X(T) is therefore a 
stationary Gaussian process with correlation function 
(X(T)X(T')) = /(T-T'). 

The simulations were performed by constructing Gaus- 
sian (pseudo-)white noise directly in the frequency do- 



main, normalizing by the appropriate y /(^), then fast- 
Fourier-transforming back to the time domain. The dis- 
tribution of intervals between zeros, and hence the per- 
sistence probability, may then be measured directly. The 
resulting process X(T) is periodic, but this is not ex- 
pected to affect the results provided the period N8T is 
sufficiently long, where N is the number of lattice sites 
used and ST is the time increment between the lattice 
sites. It is desirable for N to be as large as possible, 
consistent with computer memory limitations — typically 
N = 2 19 or 2 20 . The timestep ST must be sufficiently 
small for the short-time behaviour of the correlation func- 
tion to be correctly simulated, but also sufficiently large 
that the period of the process is not too small. Typi- 
cal values for ST were in the range 10 _4 -10 -2 . Several 
different values of ST were used for each j3, to check for 
consistency, and the results were in each case averaged 
over several thousand independent samples. 

This method works best for processes that are 
"smooth" . The d ensity o f zeros for a stationary Gaus- 
sian process is v/— /"(0) /n jLO), which is only finite if 
f(T) = 1 — 0(T 2 ). However, the correlation functions 
for the processes under consideration behave like ( |fl| ) 
at short time, so they have an infinite density of zeros. 
Any finite discretization scheme will therefore necessar- 
ily miss out on a large (!) number of zeros, and will 
presumably overestimate the persistence probability and 
hence underestimate 9 (this drawback is not present if 
Eqn. ([!]) is simulated directly, since ST = St/t effectively 
becomes zero as t — > oo). Nevertheless, the simulations 
were found to be well-behaved when (3 was greater than 
0.5, with consistent values of 9 for different values of ST. 
However, when [3 was less than 0.5 it was found to be 
increasingly difficult to observe such convergence before 
finite size effects became apparent. Convergence was not 
achieved for j3 < 0.45. Figure 4 shows the persistence 
probability for both the steady state and the transient 
processes with (3 = 0.45 and (3 = 0.75, for two values 
of ST in each case. The agreement of the data for the 
different values of ST is better for the larger value of (3. 
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FIG. 4. (a) Persistence probability obtained from simula- 
tions of the equivalent stationary Gaussian process for both 
the transient and the steady-state case, with f3 — 0.45. Two 
values of the time increment ST are shown, (b) Same as (a) 
for = 0.75. 



A summary of the measured values of 9 for different 
values of (3 is to be found in Table I. The quoted errors 
are subjective estimates based on the consistency of the 
results for different values of ST, and are smaller for the 
larger values of (3 due to the process being "smoother" . 
For (3 > 0.5, the estimated values of 9s agree well with 
the prediction 1— (3, while the result for t9 when (3 — 0.55 
agrees well with the perturbative prediction w 0.409 from 
Eqn. (H). For (3 = 0.45, the prediction 9 S = 1-/3 lies 
outside the quoted error bars, reflecting the difficulty of 
assessing the systematic errors in the simulations, while 
the result for 9o is consistent with the prediction « 0.591 
of the perturbation theory. 
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V. THE PERSISTENCE EXPONENT OF 
FRACTIONAL BROWNIAN MOTION 



length t within the interval [ti, t\ + L] is proportional to 
q(r), and can be written as 



The numerical results presented in the preceding sec- 
tion, as well as the pertubative calculation of Section 
III, clearly demonstrate the validity of the identity (|^), 
6s = 1 — fl, over the whole range < < 1. Since 
the special property of the steady state process which is 
responsible for this simple result is obscured by the map- 
ping performed in Section II, we now return to the orig- 
inal, unsealed process in time t. The crucial observation 
is that the limiting process hin^^oo H(x, t; to) defined in 
(||,[)]) has, by construction, stationary increments, in the 
sense that 



n(r) =n {L)T- {1+Ss \ 



(47) 



a 2 (t,t') = lim ((#(*,«; to) - H(x,t';t )) 2 ) 
= 2K\t-t'\ 20 



(43) 



depends only on t — t' . The power law behaviour of the 
variance (|43| ) is the defining property of the fractional 
Brownian motion (ffim) introduced by Mandelbrot and 
van Ness ||, and identifies (3 as the 'Hurst exponent' of 
the process. 

The first return statistics of fBm has been addressed 
previously in the literature ]l9|-pl"|, and analytic argu- 
ments as well as numerical simulations supporting the 
relation (^) have been presented. It seems that the re- 
lation in fact applies more broadly, to gen eral self-affinc 
processes which need not be Gaussian po| , p2[ . For com- 
pleteness we provide in the following a simple derivation 
along the lines of ||-||. 

We use H{t) as a short hand for the fBm limit of 
H(x,t;to) for to — > oo. Let Hi = H(ti), and define 
p(r) as the probability that H has returned to the level 
Hi (not necessarily for the first time) at time t\ + r. 
Obviously, using (p3[) we have 



P(r) 



1 



-0 



27rcr(r) 



(44) 



The set of level crossings becomes sparser with increas- 
ing distance from any given crossing. It is 'fractal' in the 
sense that the density, viewed from a point on the set, 
decays as r~ (1 ~ D) with D = 1 - 3. 

We now relate the decay (|44|) to that of the per- 
sistence probability ps(t)- Consider a time interval of 
length L 1. According to (|4^), the total number 
N(L) of returns to the level Hi = H(ti) in the inter- 
val ti < t < ti + L is of the order 



Now let 



N(L) ~ L 1_/3 . 



q(r) ~ -dps/dr ~ T ~ (1+( > s) 



(45) 



(46) 



denote the probability distribution of time intervals be- 
tween level crossings. The number n(r) of intervals of 



where the prefactor no(L) is fixed by the requirement 
that the total length of all intervals should equal L, i.e. 



dr to(t) ~ L. 



(48) 



This gives no(L) ~ L 6s , and since N(L) ~ n (L) com- 
parison with (45) yields the desired relation (g). 

In fijj an essentially equivalent argument was given, 
however it was also assumed that the intervals between 
crossings are independent, which is clearly not true due 
to the strongly non-Markovian character of the fBm ||. 
Here we see that what is required is not independence, 
but only stationarity of interval spacings (that is, of in- 
crements of H). The latter property does not hold for 
the transient problem, where the probability of an inter- 
val between subsequent crossings depends not only on its 
length, but also on its position on the time axis. For the 
transient process the variance of temporal increments can 
be written in a form analogous to (|43j) , 

a 2 (i, t ') = ((h(x, t) - h{x, t')) 2 ) = 2Kap(t'/t)\t - t'\ 2! \ 

(49) 



where the positive function ap interpolates monotonically 
between the limiting values ap(0) — 2 2 ' 3 " 1 and ap(l) = 1. 
This would appear to be a rather benign (scale-invariant) 
'deformation' of the fBm, however, as we have seen, the 
effects on the persistence probability are rather dramatic 
for small (3. Similar considerations may apply to the 
Riemann-Liouvillc fractional Brownian motion [p| , which 



shares the same kind of temporal inhomogeneity 27|. 

It is worth pointing out that stationary Gaussian pro- 
cesses with a short time singularity f(T) ~ 1 — C(|T| Q ), 
T — > (compare to ([[l])) have level crossing sets which 
are fractal, in the mathematical sense |23| , with Haus- 
dorff dimension D = 1 — a/2 [^^H). However, as was 
emphasized by Barbe p6f , this result only describes the 
short time structure of the process; a suitably defined 
scale-dependent dimension of the coarse-grained set al- 
ways tends to unity for large scales, since the coarse- 
grained density of crossings is finite. In other words, the 
crucial relation ([S]) holds for r — > but not for r — > oo. 
In the present context this implies that although the sta- 
tionary correlators /q and fg share the same short-time 
singularity (eq.(pl|)), for f this does not provide us with 
any information about the persistence exponent 8q. 
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FIG. 5. Summary of data for the transient exponent 
9q(/3), obtained from simulations of the solid-on-solid mod- 
els (squares), discretized Langevin equations (triangles) and 
the equivalent stationary process (diamonds). The bold line is 
the perturbation result (|28[), while the shaded area shows the 
range enclosed by the exact upper and lower bounds derived 
in Section VI. 



VI. EXACT NUMERICAL BOUNDS FOR 

In Section [0] several rigorous analytic bounds for 9q 
were obtained by comparing the correlator fa(T) to a 
function with known persistence exponent and using 
Slepian's theorem Having convinced ourselves, in 
the preceding section, that the expression (||) for 6>s is 
in fact exact in the whole interval < (3 < 1, we can now 
employ the same strategy to obtain further bounds by 
comparing fo(T,(3) to fs(T,f3). These bounds turn out 
to be very powerful because fo(T, (3) and fs{T,J3) share 
the same type of singularity at T — (see eqn.(O)). 

Let f A (T) and fs{T) be monotonically decreasing 
functions with the same class of analytical behaviour near 
T = 0. There will always exist numbers 6 m ; n and 6 max 
such that f A (b m ^T) < f B {T) < f A (b min T) is satisfied 
for all T. Since the persistence exponent for a process 
with correlation function f A (bT) is just b9 Al where 9 A is 
the persistence exponent for the process with correlation 
function f A (T), we can deduce from Slepian's theorem 
that b mul 9 A < 9b < b ma , x 9 A . We may therefore find up- 
per and lower bounds for the persistence exponent of a 
given process if we know the persistence exponent for an- 
other process in the same class. In particular, we may 
use the result @ for 0$ to obtain bounds for 9q- 

Finding by analytic means the most restrictive values 
of & min and 6 max such that fs(b max T,0) < f (T,(3) < 
fs(bminT, (3) is a formidable task. Our approach is to 
study the analytical behaviour near T = and T = oo 



to find values where the inequalities are satisfied in the 
vicinity of both limits, and then investigate numeri- 
cally whether the inequalities are satisfied away from 
these asymptotic regimes. The same leading small- X 1 
behaviour is obtained for fs(bT,f3) and fo(T, (3) when 

b — 2^~ 1 , and the same large-T behaviour is found when 
b = (1 - (5)/ (3 for ft < 1/2 and b = 1 for (3 > 1/2. We 
conclude that 



i i 



> 



max 
1 



V ' 

1 

( OW^ 1 I 

V ' 



1 



- 1 



for p < 1/2 
for (3 > 1/2 

for < 1/2 
for [3 > 1/2 



(50) 
(51) 



Surprisingly, we find in the majority of cases that the in- 
equalities ( p0| ) and (|5l]) are satisfied as equalities. Specifi- 
cally, in the cases where the corrections to the leading an- 
alytic behaviour has the appropriate sign to ensure that 
fs(bT,/3) is a bound for fo{T,(3) in both limits T ->• 
and T — > oo, then numerical investigation revealed it to 
be a bound for all T. We are therefore able in the follow- 
ing cases to quote the most restrictive bounds and their 
region of applicability in analytical form, although the 
validity of the bounds has only been established numeri- 
cally. 



> 



for </3 <(3i(= 0.1366...) (52) 

+ - 1 for ft < (3 < \ (53) 

for #2 (= 0.1549...) </3< | (54) 

for (3 > \ (55) 



% > (1-/3)2** 
(I -/?) 2 

%<l-(3 



The two critical values /?i 
0.1549 



0.1366... and (3 2 



correspond to the solutions in ]0, 0.5[ of ^- — 

1 = and j3 2 = 2 2l32 ~ 3 respectively. Note that @ 

coincides with the rigorous result in Eq. ( |l8| ) 

For other values of f3, the sign of the leading corrections 
implies that (^o|) and ( fill) is only satisfied as an inequal- 
ity, and the best value of the bound has to be obtained 
numerically by finding the value of b where fs{bT,f3) 
touches fo(T,(3) at a point. 

It is also possible to obtain bounds on 9^((3) by com- 
paring f (T,(3) with f s {bT,(3'), where /3' ^ f3. Consid- 
eration of the behaviour at small-T shows that a lower 
bound on 9q may be obtained when (3' > (3, whereas an 
upper bound may be obtained for (3 1 < (3. In the ma- 
jority of cases, it can be shown that the most restrictive 
bounds are in fact obtained by (3' — (3. For instance, con- 
sideration of the large-T behaviour for < (3 < j3\ shows 

that 6 min < i-1, so 9 min < (i - l) (1-/3). Therefore, 

since this inequality is satisfied as an equality for j3' = (3 
(see (|52|)), and /?' > /?, the best bound is obtained when 
(3' = (3. Similarly, the inequalities Q) and (E^ can be 







shown to be the best obtainable by this method. How- 
ever, for Pi < [3 < 1/2, a perturbative consideration with 
(3' = j3 + e shows that a larger value of 9 m [ n may be ob- 
tained for e small and positive, so (|53| ) is not the best 
bound that may be obtained. A numerical investigation 
shows that the optimal value of /?' is nevertheless very 
close to (3, and an improvement in 9 m in of no more than 
w 3% was so obtained. For (3 outside the ranges quoted 
in (p 7 2|-p5[), it was found by numerical investigation that 
(3 = f3' gave the best bound. 

Numerical values of these bounds are listed in Table I, 
for comparison with the simulation data. The contents 
of the last three columns of this table are also plotted in 
Figure 5. The upper and lower bounds are (perhaps) sur- 
prisingly close together. Notice that the lower bound for 
9q when j3 = 0.125 is 6.125, whereas the discrete solid- 
on-solid simulations yielded the inconclusive value > 3.3. 
All the other data are consistent with the bounds, within 
numerical error. It is interesting to note that the data, 
as well as the exact perturbation theory result, tend to 
lie much closer to the lower than the upper bound. 



fluctuations on crystal surfaces observed with the scan- 
ning tunneling microscope [p8| . From such measurements 
the autocorrelation function of the step position can be 
extracted, and different values of (3 have been observed, 
reflecting different dominant mass transport mechanisms 
p9| . Thus it seems that, perhaps with a slight refinement 
of the observation techniques, the first passage statistics 
of a fluctuating step may also be accessible to experi- 
ments. 
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VII. CONCLUSIONS 



APPENDIX: DERIVATION OF SOME 
EXPONENT INEQUALITIES 



In this paper we have investigated the first passage 
statistics for a one-parameter family of non-Markovian, 
Gaussian stochastic processes which arise in the context 
of interface motion. We have identified two persistence 
exponents describing the short time (transient) and long 
time (steady state) regimes, respectively. 

For the steady state exponent the previously conjec- 
tured relation 9s — 1 — (3 was confirmed. While this 
relation follows from simple scaling arguments applied to 
the original process, it is rather surprising when viewed 
from the perspective of the equivalent stationary process 
with autocorrelation function ([To|): It provides the exact 
decay exponent for a family of correlators whose class 
a = 2(3 (in the sense of Slepian 0) covers the whole in- 
terval a €]0, 2[; previous exact results were restricted to 
a = 1 and a — 2 pj. We have demonstrated in Section 
VI how this can be exploited to obtain accurate upper 
and lower bounds for other processes within the same 
class. 

Estimates for the nontrivial transient exponent 9q were 
obtained using a variety of analytic, exact and perturba- 
tive approaches, as well as from simulations. The numeri- 
cal techniques - direct simulation of interface models and 
construction of realizations of the equivalent stochastic 
process, respectively - are complementary, in the sense 
that the former is restricted to the regime (3 < 1/2, while 
the latter gives the most accurate results for (3 > 1/2. 
The results summarized in Figure 5 provide a rather com- 
plete picture of the function 6q{(3). 

Finally, we briefly comment on a possible experimental 
realization of our work. Langevin equations of the type 
(pj) are now widely used to describe time-dependent step 



To establish (|i7| ) we need to show that fs(T) > 
exp(-|T|/2) for all T, provided that (3 > 1/2. To this 
end we rewrite dlG) in the form 



/5(T) = > T I[1 



-2/3ITI 



-\T\\W 



(Al) 



and notice that for /3 > 1/2, [1 - exp(-|T|)] 2/3 < 1 
exp(-|T|). Thus 



fs{T)>l[e-P\ T \ + e-Wm]> 



\T\/2 



(A2) 



where the last inequality follows from the fact that the 
expression in the square brackets is an increasing func- 
tion of (3 for (3 > 1/2. 

Next we consider the relations (|l^). We express ([?]) in 
the form 



where the function g is given by 

g(y) = (l + yff } ~(l~y)^. 



(A3) 



(A4) 



Taking two derivatives with respect to y it is seen that 
g" > for (3 < 1/2 and g" < for (3 > 1/2. Since 
g(0) = and <?(1) = 2 2/3 always, it follows that g(y) 
is bounded by the linear function 2 2 @y, from above for 
f3 < 1/2 and from below for (3 > 1/2. Inserted back into 
|) this implies 



fo(T) < e - (1 ^ )|T| for (3 < 1/2 
/o(T) > e-^^l for f3 > 1/2, 



(A5) 
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and Jl8| ) follows by applying Slepian's theorem M in con- 
junction with the fact that 9 — A for purely exponential 
(Markovian) correlators. 

The inequalities ([l9]) are a little more subtle to prove. 
Let us first consider the case 1/2 < < 1. We need 
to prove that / (T) > f s (T) for all T. Then the rela- 
tion 8o < 9s will follow from Slepian's theorem. De- 
noting y — e~' T ' and using the expressions of fo(T) 
and fs(T), we then need to prove that the function 
F(y) = (1 + y) 2/3 + {a- 1)(1 - y) 2/3 - a(y 2 ^ + 1) (where 
a = 2 2 ' 3 ~ 1 ) is positive for all < y < 1. 

First we note, by simple Taylor expansion around y = 
and y = 1, that F(y) > for y close to and 1. The first 
derivative, F'(y) starts at the positive value 2/3(2 — a) at 
y = and approaches from the negative side as y — ► 1~ . 
The second derivative F"{y) starts at — oo as y — > + and 
approaches +oo as y — > l - . We first show that F"(y) is 
a monotonically increasing function of y in < y < 1. 

To establish this, we consider the third derivative, 
F"'(y) = 2/3(2/3 - l)(2/3 - 2)G(y), where 



= (1 + j/) 



2/3-3 



(a-l)(l-») 



2/3-3 



ay 



2/3-3 



(A6) 



Now, since (1 + y)( 2 ^" 3 ) < yC 2 ?- 3 ) for all < y < 1, 
we have G(y) < — (a — l)[y 2/5 ~ 3 + (1 — y) 2/3 ~ 3 ] implying 
G(y) < for < y < 1. Since 1/2 < < 1, it follows 
that F"'(y) > for all < y < 1. Therefore, F"(y) is 
a monotonically increasing function of y for < y < 1 
and hence crosses zero only once in the interval [0,1]. 
This implies that the first derivative F'{y) has one sin- 
gle extremum in [0, 1]. However, since F'(y) starts from 
a positive value at y = and approaches from the 
negative side as y — ► 1~, this single extremum must be 
a minimum. Therefore, F'(y) crosses zero only once in 
[0, 1] implying that the function F(y) has only a single 
extremum in [0, 1]. Since, F(y) for y — ► + and y — > 1~ , 
this must be a maximum. Furthermore, -F(y) can not 
cross zero in [0, 1] because that would imply more than 
one extremum which is ruled out. This therefore proves 
that /o(T) > f s (T) for all T for > 1/2 and hence, 
using Slepian's theorem, 9q < 9s- Using similar argu- 
ments, it is easy to see that the reverse, 9q > #S is true 
for < 1/2. 



Finally we prove the inequality (14) which relates the 
values of 9q for two different exponents and 0' > 0, 
subject to an additional constraint to be specified be- 
low. In the same spirit as above, one can show that after 
defining 7 = (1 — /3)/(l — 0') > 1, one obtains 



fo(T,0)<fo(lT,0') 

for 0<0' and 02~ 2f) < fiT 2 ?' . 



(A7) 



Both functions in (A7) decay exponentially at large time 
with the same decay rate Ao = 1 — = 7(1 — 0'), such 
that their ratio approaches a constant when T — > 00. 
The last condition 02~ 2/3 < /3'2" 2 ' 3 ' expresses that the 
limit of this ratio must be less than unity. As fo{T, 0) < 



foi'jT, 0') in the vicinity of T — (this is just a con- 
sequence of < 0'), and as a careful study shows that 
-^p [fo (T,0)/ f (7T, 0')] has at most one zero in the range 
]0,oo[, we conclude that f (T, 0)/ 'fo^T, 0') < 1 if and 
only if the limit of this ratio for T — > is less than 
unity. In practice, the last condition in (A7) express- 
ing this constraint can be violated only if > 1/2 and 
0' > (21H2)- 1 w 0.7213475.... 

Using Slepian's theorem, and the fact that the per- 
sistence exponent associated with /o(7T, 0) is exactly 7 
times the persistence exponent associated with fo (T, 0) 
we arrive at the inequality (|lj) which holds under 
the conditions stated in (Wfj ). Setting or 0' equal to 
1/2 (keeping < 0'), eq.(|14|) reduces to the bounds jig)). 

For 0' > > 1/2 the inequality ( |l4| ) comes rather close 
to being satisfied as an equality. For example, setting = 
0.55, 0' = 0.85 eq.dH) requires that 9 O (0)/9 O (0') > 3, 
while the numerical data in Table I yield 9q(0) / 9o(0') = 
3.39. The only pair of values in Table I which violates the 
inequality Q is {0,0') = (0.75,0.85). Since for these 
values the condition /32~ 2/3 < 0'2~ 2 ^ is also violated, 
this may be taken as an indication that the numerical 
estimates for > 1/2 are rather accurate. 
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p 


0s 


1-/3 


00 




^max 


0.125 ( * ) 


0.86 ±0.02 


0.875 


> 3.3 


6.125 


7.359 


0.2 


0.788 ±0.01 


0.8 


2.0 ±0.1 


2.333 


3.200 


0.25W 
0.25 


0.754 ±0.01 
0.740 ±0.01 


0.75 
0.75 


1.6 ±0.15 
1.55 ±0.02 


1.547 
1.547 


2.250 
2.250 


0.3 


0.69 ±0.01 


0.7 


1.10 ±0.05 


1.141 


1.633 


0.375^ 
0.375 


0.635 ±0.01 
0.625 ±0.01 


0.625 
0.625 


0.84 ±0.01 
0.85 ±0.01 


0.801 
0.801 


1.042 
1.042 


0.4 


0.60 ±0.01 


0.6 


0.76 ±0.1 


0.723 


0.900 


0.45 


0.53 ±0.01 


0.55 


0.58 ±0.02 


0.598 


0.672 


0.55 


0.44 ±0.01 


0.45 


0.41 ±0.01 


0.415 


0.450 


0.6 


0.39 ±0.01 


0.4 


0.35 ±0.01 


0.348 


0.400 


0.65 


0.346 ±0.005 


0.35 


0.295 ±0.005 


0.289 


0.350 


0.75 


0.247 ±0.005 


0.25 


0.201 ±0.005 


0.191 


0.250 


0.85 


0.150 ±0.005 


0.15 


0.121 ±0.005 


0.107 


0.150 



Table I. Numerical estimates of the persistence exponents 63 and 6q, compared to the fractional Brownian motion 
result 9s = 1 — /3 (third column) and to the optimal bounds # m in and 9 max derived in Section VI (last two columns) . 
With the exception of the values marked with an asterisk'*), which were obtained using discrete solid-on-solid models 
(Section IV.A), the data for (3 < 0.40 are taken from simulations of discretized Langevin equations (Section IV. B) 
while those for /3 > 0.45 were generated using the equivalent stationary Gaussian process (Section IV. C). In all cases 
the error bars reflect a subjective estimate of systematic errors. 
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